\chapter{Introduction}

\section{Image Segmentation}
Image segmentation is the task of splitting a digital image into one or more regions of interest. It is a fundamental problem in computer vision and many different methods, each with their own advantages and disadvantages, exist for the task. Image segmentation is a particularly difficult task for several reasons. Firstly, the ambiguous nature of splitting up images into objects of interest requires a trade off between making algorithms more generalized or having many user specified parameters. Secondly, image artifacts such as noise, inhomogeneity, acquisition artifacts and poor contrast, are very difficult to account for in segmentation algorithms without a high level of interactivity from the user. 

In this report, segmentation is discussed in a medical imaging context however the proposed algorithm could equally be used in general purpose segmentations. Segmented images are typically used as the input for applications such as classification, shape analysis and measurement. In medical image processing, segmented images are used for studying anatomical structures, diagnosis and assisting in surgical planning. Before computational algorithms existed, segmentation of medical images was a tedious process, performed by hand by clinical experts. This was a fairly accurate, yet slow process. These expert segmentations form the gold standard with which to validate algorithmic segmentations.

As there is no single general solution to image segmentation problems, several techniques exist each of which has its own strengths and limitations. Some common techniques are thresholding based segmentation (discussed in Section \ref{thresholding}), region based segmentation (i.e. \textit{region growing}), edge based segmentation and deformable active contour models (such the \textit{snakes} \cite{kass1988sac} and geodesic active contour models).

\section{Motivation}
In biomedical image segmentation, there is great interest segmenting three dimensional anatomical structures from CT, MRI and PET data. General purpose segmentation on the other hand is typically a two dimensional problem. As a result the data sets used to segment biomedical images are much larger ($O(n^3)$) than those used otherwise ($O(n^2)$), resulting in computation times an order of magnitude greater (although this is not true in all cases). In addition to the time complexity there is also significant storage complexity, for example a very large MRI volume could have as many as $512^3$ elements or more, which in combination with large element sizes may result in memory errors on current available hardware when seeking a large area of contiguous memory to write this to. 

Furthermore, although the level set approach used in this project is very powerful (its advantages are discussed in Section \ref{levelsetmethod}), the principal disadvantage of using level sets is that they are relatively slow to compute.

As a result, processing times in a clinical setting may be impractical with more computationally demanding segmentations requiring many minutes or even hours. This provides us with the motivation to accelerate these segmentations by using a parallel algorithm executed on the GPU (Graphics Processing Unit).


\section{Parallel Processing}
The algorithms for processing level sets have vast parallelization potential. Section \ref{levelsetalgorithm} details the algorithms used to discretize the level set equation and Section \ref{parallel} discusses how these can be executed on graphics hardware.
	\subsection{GPGPU}
General purpose computation on graphics processing units (GPGPU) is the technique of using graphics hardware to compute applications typically handled by the central processing unit (CPU). Graphics cards over the past two decades have been required to become highly efficient at rendering increasingly complex 3D scenes at high frame rates. This has forced their architecture to be massively parallel in order to compute graphics faster than general purpose CPUs.

Compared to a CPU, a GPU features many more transistors on the control path due to the lower number of control instructions required. Memory is optimized for throughput and not latency, with strict access patterns. GPUs are not optimized for general purpose programs, and they do not feature the complex instruction sets, or branch control of the modern CPU. Finally, although current high-performance CPUs feature multiple cores for limited parallel processing, GPUs are arguably a more attractive option in terms of lower price and power usage.

The advent of GPGPU programming came with programmable shader units that allowed the programmer to write small programs at each pixel or vertex in the rendering pipeline. In order to program the shader units, shading languages had to be used in conjunction with graphics APIs such as DirectX and OpenGL. NVIDIA developed the high-level shading language Cg to assist in programming shaders, however it still required knowledge of graphics APIs. More recently, languages have been developed that allow the programmer to implement algorithms without any knowledge of graphics APIs or architectures. One such language is NVIDIA CUDA, and is the language chosen for the optimizations in this project.

	\subsection{CUDA}
Compute Unified Device Architecture, or CUDA, is NVIDIA's GPGPU technology that allows for programming of the GPU without any graphics knowledge. The C language model has at its core three key abstractions, from \cite{cuda}: a hierarchy of thread groups (to allow for transparent scalability), shared memories (allowing access to low-latency cached memory), and barrier synchronization (to prevent \textit{race conditions}). This breaks the task of parallelization into three sub problems, which allows for language expressivity when threads cooperate, and scalability when extended to multiple processor cores.

		\subsubsection{Framework}
CUDA uses extends C by allowing a programmer to write \textit{kernels} that when invoked execute thousands of lightweight identical threads in parallel. CUDA arranges these threads into a hierarchy of blocks and grids, as can be seen in Figure \ref{fig:cudathreads} allowing for runtime transparent scaling of code to different GPUs. The threads are identified by their location within the grid and block, making CUDA perfectly suited for tasks such as image processing where each thread is easily assigned to an individual \textit{pixel} (picture element) or \textit{voxel} (volume element).

\begin{figure}[p]
	\centering
		\includegraphics[scale=0.4]{images/cudathreads.PNG}
		\caption{A grid of thread blocks. This figure is taken from \cite{cuda}}
	\label{fig:cudathreads}
\end{figure}

When writing and optimizing complex parallel code in CUDA it is often found that threads may need to cooperate, typically by sharing data or temporarily halting execution. The memory hierarchy of CUDA threads is shown in Figure \ref{fig:cudamemory}. Here it can be seen that each thread has access to: a per-thread private local memory, a per-block on-chip shared memory to share data between threads, and finally an off-chip global memory accessible to all threads within all blocks. There are also constant and texture memory spaces accessible to all threads, however these are not featured in our algorithm and so will not be discussed in any further detail.

\begin{figure}[p]
	\centering
		\includegraphics[scale=0.4]{images/cudamemory.PNG}
		\caption{The memory hierarchy of CUDA threads and blocks. This figure is taken from \cite{cuda}}
	\label{fig:cudamemory}
\end{figure}

		\subsubsection{Performance Guidelines}
There are many techniques to optimize a parallel algorithm. Firstly, the optimum block and grid sizes should be used to ensure maximum `occupancy'. Occupancy is the ratio of the number of active warps (32 parallel threads) to the maximum number of active warps supported by the GPU multiprocessor. To maximize efficiency, there is a trade off between making the occupancy high enough to ensure no multiprocessor is ever idle, and making it low enough to ensure no \textit{bank conflicts} (when two threads attempting to access the same location in shared memory have their accesses serialized).

Secondly, one of the best ways in which to optimize the parallelization is through efficient memory usage. The global memory space is not cached and therefore has a much higher latency and lower bandwidth than on-chip shared memory. Therefore it is the aim of the programmer to minimize global memory accesses. From \cite{cuda}, it recommended that each thread in a block firstly loads data from global memory to shared memory, synchronizes with all other threads within the thread block to ensure shared memory locations have been written to, processes the data, synchronizes again to ensure shared memory has been fully updated with results, and finally writes the results back to global memory coalesced.

\textit{Coalescence} is an important concept in memory management as it can speed up memory reads and stores significantly. \cite{cuda} lists the following three conditions for coalescing: ``threads must access either 32-bit words, 64-bit words, or 128-bit words'', ``all 16 words must lie in the same segment of size equal to the memory transaction size'' and ``threads must access the words in sequence''. 

%Finally, some of the strict requirements and limitations of CUDA are disappearing as devices of higher \textit{compute capability} feature incremental improvements to the core architecture. For example, devices of compute capability 1.3 have more relaxed requirements for coalescence and support for double-precision floating point accuracy not found in devices of compute capability 1.1.